宏基因组教程Metagenomics Tutorial (HUMAnN2)

Original: metagenome 宏基因组 2019-07-15

之前我们介绍了microbiome_helper提供众多宏基因组学相关的实验、分析教程。

微生物组助手——最易学的扩增子、宏基因组分析流程

今天我们先详细讲解基于HUMAnN2的有参宏基因组分析流程 Metagenomics Tutorial (HUMAnN2) https://github.com/LangilleLab/microbiome_helper/wiki/Metagenomics-Tutorial-(Humann2))

此流程在我们之前学习的相关软件基础上，进一步完善了分析流程，包括数据质控、去宿主、

中间的物种和功能组成量我们之前有过介绍，

同时还提供了各步骤间的格式转换，以及下游的统计绘图。

说明：如果使用官方的虚拟机，可完全不用安装软件和修改原版代码即可运行测试数据。本文是基于我个人系统安装了最新版本软件，因此软件名称、安装位置、数据库位置都有变化，请用户根据自身情况选择学习。

分析流程

下载测试数据

mkdir -p ~/test/mh_meta17/v2
cd ~/test/mh_meta17/v2
wget http://kronos.pharmacology.dal.ca/public_files/tutorial_datasets/mgs_tutorial_Oct2017.zip
unzip mgs_tutorial_Oct2017.zip   
cd mgs_tutorial_Oct2017
gunzip raw_data/*

了解输入文件

# 实验设计
less map.txt
# 样品数
wc -l map.txt
# 测序数据
head -n 8 raw_data/p136C_R1.fastq

输入文件1：实验设计

主要包括样品名，分组(如正常/癌症)，相关属性

Sample Id    Sample Type    Sex    Individual
p136C    Cancer    M    p136
p136N    Normal    M    p136
p143C    Cancer    M    p143

输入文件2：测序文件，每个样品1个或1对fastq/fq文件

@SRR3586062.883556
CTTGGGGCTGCTGAGCTTCATGCTCCCCTCCTGCCTCAAGGACAATAAGGAGATCTTCGACAAGCCTGCAGCAGCTCGCATCGACGCCCTCATCGCTGAGG
+
CCCFFFFFHHHHHIJJJJJJJIJIJJJJGIJDGIJEIIJIJJJJJJJJIJJJJIJJIJJJJJHHHFFFFECEEEDDDDD?BDDDDDDBDDDDDDDDBBBDD

软件安装和环境变量

export PATH=$PATH:/conda2/bin

详细安装方法见文末软件安装部分。为什么我不把安装目录永久添加在环境变量，因为conda也会影响环境变量，导致我其它工具依赖关系出错，因此特殊流程可以临时添加更安全，自己清楚在每个流程在哪里，即使添加也要放在$PATH后面才安全。

不可能存在一个环境满足所有的依赖关系，所以近几年conda, docker会非常流行。推荐安装首选conda，如果仍搞不定，再用docker(小提示，每个bioconda软件都提供docker下载，只是使用更麻烦但成功率更高)。

序列质控和去宿主

# 新版本脚本名称变更，指定bowtie2数据库、trimmomatic位置
kneaddata -h # 显示帮助
kneaddata -i raw_data_example/p144C_R1.fastq -i raw_data_example/p144C_R2.fastq -o kneaddata_out -v \
-db ~/ref/host/Homo_sapiens  \
--trimmomatic /mnt/bai/public/bin/Trimmomatic-0.36/ --trimmomatic-options "SLIDINGWINDOW:4:20 MINLEN:50" \
-t 9 --bowtie2-options "--very-sensitive --dovetail" --remove-intermediate-output

-h 显示帮助
-v 显示运行信息，方便观察运行进展和出错找原因
-i 为输入fastq文件，双端需输两次
-o 输出结果目录
-db 指定bowtie2索引
—trimmomatic 指定质控程序位置，默认为/usr/bin/trimmomatic-0.36.jar
—trimmomatic-options 质控选项，4碱基滑窗，质量大于20，最小长度50
-t 设置线程数，不要超过9
—bowtie2-options 比对选项
—remove-intermediate-output 清理中间文件

批处理质控

你不可能实验只一个文件，一般最小2个实验组，每组三次重复，需要批量处理

这里使用parallel命令，依赖map文件中样品名，或文件列表对测序结果并行处理。注意不同版本下参数可能不同，如处理多参数时，原文虚拟机20170322版本为—links参数，而我的系统中20141022版本为—xapply参数

# 批量读取成对文件作为输入
parallel -j 3 --xapply 'echo {1} {2}'  ::: raw_data/*_R1.fastq ::: raw_data/*_R2.fastq

# 批处理样品
parallel -j 3 --xapply 'kneaddata -i {1} -i {2} -o kneaddata_out -v \
-db ~/ref/host/Homo_sapiens  \
--trimmomatic /mnt/bai/public/bin/Trimmomatic-0.36/ --trimmomatic-options "SLIDINGWINDOW:4:20 MINLEN:50" \
-t 9 --bowtie2-options "--very-sensitive --dovetail" --remove-intermediate-output' \
 ::: raw_data/*_R1.fastq ::: raw_data/*_R2.fastq

质控后结果统计

kneaddata_read_count_table --input kneaddata_out --output kneaddata_read_counts.txt

合并双端

# 清理宿主污染至指定目录备用
mkdir kneaddata_out/contam_seq
mv kneaddata_out/*_contam*fastq kneaddata_out/contam_seq

# 本质上只合并了1/2
concat_paired_end.pl -p 9 --no_R_match -o cat_reads kneaddata_out/*_paired_*.fastq

# 可选方法2：我感觉是合并4个文件，应该重命名序列ID，不则双端序列重名字，shell命令合并单样品
for i in `tail -n+2 map.txt | cut -f 1`;do \
cat kneaddata_out/${i}_R1_kneaddata_paired_1.fastq kneaddata_out/${i}_R1_kneaddata_paired_2.fastq kneaddata_out/${i}_R1_kneaddata_unmatched_1.fastq kneaddata_out/${i}_R1_kneaddata_unmatched_2.fastq | awk '{if(NR%4==1) print "@"NR; else print $0}' > cat_reads/${i}.fastq; done

计算功能和代谢通路

humman2是一套分析流程，它包括调用metaphlan流程来分析物种组成，和自身分析功能基因和代谢通路组成。安装方法见本文软件安装部分。

humann2 --threads 9 --input cat_reads/p144C.fastq --output humann2_out/

多样品批量并行处理，-j参数请尽量小于你你电脑核心数(线程数的一半)，否则效率反而会低且系统卡顿。

parallel -j 9 'humann2 --threads 9 --input {} --output humann2_out/{/.}' ::: cat_reads/*fastq

多样品物种和功能组成合并为矩阵/表

merge_metaphlan_tables.py precalculated/metaphlan2_out/*tsv > metaphlan2_merged.tsv

# 转换为spf格式方便stamp分析
metaphlan_to_stamp.pl metaphlan2_merged.tsv > metaphlan2_merged.spf
# 删除样品名中多余字符，和株水平行
sed -i 's/_metaphlan_bugs_list//g' metaphlan2_merged.spf 
# 删除株行：新数据库新增行
sed -i '/t__[^\t]*\t/d' metaphlan2_merged.spf

STAMP软件统计绘图

可以对上述结果使用STAMP打开，鼠标点点可完成常见统计分析的可视化。主页 http://kiwi.cs.dal.ca/Software/STAMP

可参考我们之前的教程

微生物组间差异分析神器-STAMP

可进一步在软件主页学习第一手英文教程。

整理humann2功能结果

humann2_join_tables --input precalculated/humann2_out/ --file_name pathabundance --output humann2_pathabundance.tsv

# 标准化为百分比
humann2_renorm_table --input humann2_pathabundance.tsv --units relab --output humann2_pathabundance_relab.tsv

# 分层结果
humann2_split_stratified_table --input humann2_pathabundance_relab.tsv --output ./

# 筛选某个通路结果
head -n 1 humann2_pathabundance_relab_unstratified.tsv >> humann2_pathabundance_relab_LACTOSECAT-PWY.tsv
grep "LACTOSECAT-PWY" humann2_pathabundance_relab_unstratified.tsv >> humann2_pathabundance_relab_LACTOSECAT-PWY.tsv

# 格式化为stamp输入
sed 's/_Abundance//g' humann2_pathabundance_relab_LACTOSECAT-PWY.tsv >  humann2_pathabundance_relab_LACTOSECAT-PWY.spf
sed -i 's/# Pathway/MetaCyc_pathway/' humann2_pathabundance_relab_LACTOSECAT-PWY.spf

使用STAMP分析结果

软件安装

本流程虽然只有简单的几步，但是背后却是几十步，需要的软件极多，也许安装和配置的复杂程序，反倒是分析流程难度的几倍。

作者提供了Virtalbox虚拟机 https://github.com/LangilleLab/microbiome_helper/wiki/Microbiome-Helper-Virtual-Box，有20GB大小，可以直接学习，但灵活性不高，运行效率也低。有服务器的朋友应该更喜欢正常安装软件使用更高效方便。

依赖软件清单如下，这里仅对我使用的系统没有的软件进行安装说明，其它软件请读者自行查看软件官方帮助安装。

https://github.com/LangilleLab/microbiome_helper/wiki/Requirements

流程相关脚本

为完成流程各软件间的衔接，需要写很多脚本，作者整理了几十个常用脚本工具，可下载并添加环境

cd ~/github
git clone git@github.com:LangilleLab/microbiome_helper.git
echo 'export PATH=$PATH:~/github/microbiome_helper' >> ~/.bashrc

# 安装perl模块
cpan
install File::Basename Getopt::Long List::Util Parallel::ForkManager Pod::Usage

质控和去宿主kneadData

https://bitbucket.org/biobakery/kneaddata/wiki/Home

KneadData是一款宏基因组和宏转录组测序数据质控的流程，其主要功能包括使用Trimmomatic序列质控，bowtie2比对至宿主和PhiX基因组去除宿主和测序平衡序列。

#以下三个方式均不可用
#pip install kneaddata # 无权限 
#pip install kneaddata --user # 不满足依赖关系
#pip3 install kneaddata --user # 安装成功，位于~/.local/lib/python3.5/site-packages/kneaddata，配置了依赖软件和数据库，在bowtie步仍报错

# /conda2安装
/conda2/bin/conda install kneaddata

# 列出主要命令
ll /conda2/bin/kneaddata*

echo 'export PATH=$PATH:/conda2/bin' >> ~/.bashrc

# 或者临时添加conda2加入环境变量
export PATH=$PATH:/conda2/bin

# 如果你还不可用，使用docker吧，https://quay.io/repository/biocontainers/kneaddata，每个conda都有docker可下载，但也需要更大权限

查看可用数据

kneaddata_database

下载人类和小鼠基因组数据库

mkdir -p ~/ref/host
cd ~/ref/host
kneaddata_database --download human_genome bowtie2 ./
kneaddata_database --download mouse_C57BL bowtie2 ./

创建定制数据库

本质上就是把指定的基因组建一个bowtie2的索引

https://bitbucket.org/biobakery/kneaddata/wiki/Home#markdown-header-create-a-custom-database

以下载的EMBL plants上的拟南芥基因组为例

# bowtie2-build <reference> <db-name>
cd ~/ref/ath
bowtie2-build Arabidopsis_thaliana.TAIR10.dna.toplevel.fa bt2ebi

humann2安装

详见之前的文章

HUMAnN2：人类微生物组统一代谢网络分析2

系列教程：微生物组入门 Biostar 微生物组宏基因组

专业技能：生信宝典学术图表高分文章不可或缺的人

一文读懂：宏基因组寄生虫益处进化树

必备技能：提问搜索 Endnote

文献阅读热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读分析流程统计绘图

16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun

在线工具：16S预测培养基生信绘图

科研经验：云笔记云协作公众号

编程模板 Shell R Perl

生物科普生命大跃进细胞暗战人体奥秘

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外150+ PI，1500+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读

”FAN某”的离婚财产分割判决书（全文）

”FAN某”的离婚财产分割判决书（全文）

公益慈善｜“翼行天下一生守护”慈善项目捐赠仪式圆满举行！

哈里斯女粉搞4B运动、毒杀丈夫，回旋镖能否让美国“血流成河”

比国产光刻机更重要的IPO要来了！

宏基因组教程Metagenomics Tutorial (HUMAnN2)

分析流程

下载测试数据

了解输入文件

软件安装和环境变量

序列质控和去宿主

质控后结果统计

合并双端

计算功能和代谢通路

多样品物种和功能组成合并为矩阵/表

STAMP软件统计绘图

整理humann2功能结果

软件安装

流程相关脚本

质控和去宿主kneadData

humann2安装

猜你喜欢

写在后面

您可能也对以下帖子感兴趣

”FAN某”的离婚财产分割判决书（全文）

”FAN某”的离婚财产分割判决书（全文）

公益慈善｜“翼行天下 一生守护”慈善项目捐赠仪式圆满举行！

哈里斯女粉搞4B运动、毒杀丈夫，回旋镖能否让美国“血流成河”

比国产光刻机更重要的IPO要来了！

生成图片，分享到微信朋友圈

宏基因组教程Metagenomics Tutorial (HUMAnN2)

分析流程

下载测试数据

了解输入文件

软件安装和环境变量

序列质控和去宿主

质控后结果统计

合并双端

计算功能和代谢通路

多样品物种和功能组成合并为矩阵/表

STAMP软件统计绘图

整理humann2功能结果

软件安装

流程相关脚本

质控和去宿主kneadData

humann2安装

猜你喜欢

写在后面

您可能也对以下帖子感兴趣

公益慈善｜“翼行天下一生守护”慈善项目捐赠仪式圆满举行！